Evolution and consequences of individual responses during the COVID-19 outbreak

In a long-lasting major disease outbreak such as that of COVID-19, the challenge for public health authorities is to keep people motivated and keen on following safety guidelines. In this study, a compartmental model with a heterogeneous transmission rate (based on awareness) is utilized to hypothesize about the public adoption of preventive guidelines. Three subsequent outbreaks in South Korea, Pakistan, and Japan were analyzed as case studies. The transmission, behavior change, and behavioral change ease rates of the disease were measured in these countries. The parameters were estimated using the maximum likelihood method with an additional identifiability analysis performed to determine the uniqueness of the estimated parameters for quantitatively comparing them during the first three waves of COVID-19. The mathematical analysis and simulation results show that individual responses had a significant effect on the outbreak. Individuals declining to follow the public health guidelines in Korea and Japan between the second and third waves contributed to making the third peak the highest of the three peaks. In Pakistan, however, individual responses to following public health guidelines were maintained between the second and third waves, resulting in the third peak being lower than the first, rather than being associated with the highest transmission rate. Thus, maintaining a high level of awareness is critical for containing the spread. Improvised public health campaigns are recommended to sustain individual attention and maintain a high level of awareness.


Introduction
The coronavirus disease  is an infectious disease caused by the SARS-CoV-2 virus. Most people infected with the virus experience mild-to-moderate respiratory symptoms and recover without special treatment. However, some patients can become seriously ill and may require medical attention. Older people and those with underlying medical conditions such as cardiovascular diseases, diabetes, chronic respiratory illnesses, or cancer, are more likely to develop serious illnesses following COVID-19 infection. Anyone can get infected with COVID-19 and serious consequent health issues, including death, can never be ruled out with certainty. Chinese authorities have identified this outbreak as a new coronavirus, differing completely from other coronaviruses previously encountered among humans [1]. The novel a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 We extend the behavior-disease modeling approach proposed in [28]. The model was fitted to each of the first three COVID-19 waves that occurred in South Korea, Pakistan, and Japan. It is important to accurately associate unknown parameters with the values that relate the model to the data. For a comprehensive analysis of the modeling results, it is also necessary to determine the reliability of the parameter estimates. The identifiability analysis guarantees that the estimated parameters are uniquely determined. If the parameters are identifiable, the parameter estimation is independent of the initial guess of the unknown parameters. Ensuring the identifiability of the parameters, a model was fitted to explore the change in individual responses to public health guidelines, demonstrating that the tendency to adopt preventive measures inspired by public health initiatives has decreased over time in Korea and Japan but is maintained in Pakistan.

Mathematical modeling
As a modification to the model proposed in [28], heterogeneity in transmission is assumed based on the level of awareness of individuals. The proposed model consists of five classes, where the total population N(t) was partitioned into classes S(t), S f (t), I(t), Q(t), and R(t). These classes denote susceptible, behavior-changed susceptible, infectious, active (quarantined or hospitalized), and recovered individuals, respectively. A schematic representation of this model is shown in Fig 1. The system of equations associated with the aforementioned diagram can be expressed as: In the susceptible class S, individuals become infected with the transmission rate β and move to the infectious class I. When the number of infected cases increases in society, fear of the disease prevails in susceptible individuals. With the behavior change rate (transmission rate of fear) β f , susceptible individuals move to the behavior-changed susceptible class S f . Because of this behavioral change, the transmission rate in S f is lowered by a fraction r β (0�r β �1). Therefore, individuals from S f move to infectious class I with a lower transmission rate, r β β than those from S. Strictly following social distancing, mask-wearing, and avoiding social gatherings are considered as behavioral changes in the proposed model. After confirmation, infectious individuals move to quarantine class Q with a confirmation rate γ and confirmed infected individuals Q, move to the recovered class with a rate of recovery α. Individuals in the model are recruited to susceptible class S with a constant migration rate Λ, and they are removed from all classes upon natural death rate μ. In quarantined classes, individuals are also removed with disease-related death rate μ 0 . Individuals in S f influenced by the susceptible and recovered individuals, become hesitant to strictly follow the preventive measures and revert to the susceptible class at a behavioral change ease rate of μ f .

Mathematical analysis
The proposed mathematical model was analyzed to gain an understanding of the dynamics of models (1)-(5) from an epidemiological perspective.
The model exhibits three equilibria: trivial, (0,0,0,0,0) disease-free (S 0 ,0,0,0,0), and endemic equilibrium ðS � ; S � f ; I � ; Q � ; R � Þ. However, bSI N and r b bS f I N are not differentiable at the trivial equilibrium. Therefore, stability analysis at this equilibrium was not feasible using standard linearization methods. The existence and uniqueness of the remaining equilibria are described in the following propositions: Proposition 1: The system of Eqs (1)-(5) has a unique disease-free equilibrium E 0 , defined as E 0 = (S 0 ,0,0,0,0). Here, S 0 ¼ L m . Proof. For the disease-free equilibrium, the left-hand side of the system of equations was set as equal to 0, and it was assumed that there was no disease in the beginning.
Therefore, I = 0, implies Q = 0 and R = 0. From the first and second equations of the system, it is determined that: This contradicts the fact that S, N�0; therefore, S f will always be 0 in a disease-free equilibrium. This is biologically meaningful, as it is natural to assume that all susceptible individuals are fear-free in the beginning. In addition, after time t, if the disease wanes, the fear of the disease subsides in society and behavioral changes return to normal in the susceptible population.
Letting S f = 0, we obtain Eq (6) as S ¼ L m . Thus, there is a unique disease-free equilibrium

Proposition 2:
The system of Eqs (1)-(5) has a unique endemic equilibrium, E � , defined as Proof. For the endemic equilibrium, it is assumed that I6 ¼0. Using this assumption, the uniqueness of the endemic equilibrium can be easily proved.∎

Stability analysis
The basic reproduction number, R 0 , plays a crucial role in the analysis of single-population outbreaks, which are the total expected secondary cases that the primary infected individual generates. To determine R 0 , the approach defined in [29] is used. Thus, the next generation matrix, J, in our model is The eigenvalues of J are b mþg and 0, which implies that the spectral radius of J is b mþg . Hence, Lemma 1: Disease-free equilibrium E 0 is locally asymptotically stable in domain R 5 þ if and only if the basic reproduction number R 0 , is less than 1.
Proof. If there is no infectious individual, then Proposition 1 ensures the uniqueness of disease-free equilibrium E 0 . We show that this unique disease-free equilibrium E 0 , is locally asymptotically stable. The Jacobian matrix of the system of Eqs (1)-(5) at E 0 is The eigenvalues of B 1 are −μ and −(μ+μ f ) and the eigenvalues of B 4 are −μ, −(α+μ+μ 0 ) and β−(μ+γ). Therefore, if R 0 ¼ b mþg < 1, all eigenvalues of B have negative real parts. Therefore, according to the Routh-Hurwitz stability criterion, E 0 is locally asymptotically stable in R 5 þ . Lemma 2: Endemic equilibrium E � is locally asymptotically stable in domain R 5 þ if and only if the basic reproduction number R 0 is greater than 1.
Proof. We prove this result by using the center manifold theorem. Let ϕ be the bifurcation parameter. If y 1 = S, y 2 = S f , y 3 = I, y 4 = Q, y 5 = R and y = y 1 +y 2 +y 3 +y 4 +y 5 , the system of equations can be expressed as follows: The linearization matrix, D y , of the aforementioned system of equations around E 0 , when ϕ = ϕ � , is given by, It is observed that 0 is a simple eigenvalue of this matrix. Therefore, let Thus, the right eigenvector, H, is given by, Similarly, the left eigenvector of matrix D y is Therefore, coefficients a and b are defined as By substituting the values, where b = 1, when ϕ>0, a<0, and b>0. Therefore, based on Theorem 4.1 in [30], the endemic equilibrium E � is locally asymptotically stable if and only if ϕ>0. This implies that E � is locally asymptotically stable if and only if R 0 >1. ∎

Parameterization
The descriptions and values of the parameters used in the system of equations are presented in Tables 1 and 2. These parameters are used to solve the system of Eqs (1)-(5). In Table 1, the parameters are adopted from references or estimated from the data. According to World Bank, the crude birth rate per year in the Republic of Korea is 6 1000 (see [31]) and the total population is approximately 52 million. Therefore, the number of new births per year is 312,000 and the immigration rate per day is 312000 365 � 855. Similarly, the migration rate for Pakistan and Japan is 16953 and 2413, respectively.
The average life expectancy in South Korea is approximately 83 years [32]; therefore, the natural death rate in Korea is 1 83 � 1 365 per day. The average life expectancy in Pakistan is 67 years [32]. Therefore, the natural death rate in Pakistan is 1 67 � 1

365
. Similarly, the natural death rate in Japan is 1 85 � 1 365 [32]. We estimate Pakistan's disease-related death rate on average is, Deaths due to disease Confirmed cases À � ¼ 0:0202 [33]. Furthermore, the disease-related death rate in Japan is 0.0238.
In Table 2, the remaining parameters are estimated by using the maximum likelihood method. In particular, three parameters are estimated for each of the three waves. The period The timeline of the three peaks for Japan was described as follows: first wave (March 10, 2020 to June 28, 2020), second wave (June 28, 2020 to October 21, 2020), and third wave (October 21, 2020 to February 14, 2021) based on [33] (please refer to S3 Fig and S3 Table in S1 File).
For all three waves, the initial value of the infectious population is estimated using the maximum-likelihood method. Fear at the start of the disease is negligible; therefore, S f (0) = 0 in the first wave. S f is also modified to N−S−I−Q−R in our system of equations.
It is observed that with an increase in the number of unknown parameters, identifiability issues are raised. To address these issues, a code in MATLAB is executed for the first wave up to the initial value of the second wave and considering S(178) = 34.

Fitting the temporal data
In the formulation, β is the transmission rate; β f is the behavior change rate; μ f is the behavioral change ease rate; and I(0) is the initial value of the infectious population. As data on the infectious population were unavailable, estimating these unknowns directly from the data was impossible. Thus, the parameters β, β f , μ f , and I(0) were estimated using the maximum likelihood method. Time-series data were employed for active infected cases ( [33]). The details of the time-series data are provided for those countries in S1-S3 Tables in S1 File.
Let X ¼ ðx 1 ; x 2 ; x 3 ; x 4 ; x 4 Þ ¼ ðS; S f ; I; Q; RÞ 2 R 5 þ be the vector of state variables; F is the right side of the system ((1)-(5)); P = (β, β f , μ f , I(0)) is the vector of the unknowns to be estimated; Q(t, P) be the vector of observables, and Q 0 (t, P) is the observed data at t. Here, t is the time interval for the first wave (the assumptions were similar for the second and third waves). It is assumed that all Q 0 (t, P) are independent and drawn from the Poisson distribution with a mean equal to Q(t, P); then, the Poisson maximum likelihood function is Therefore, the negative log-likelihood function (NLF) reduces to Because the last term in the aforementioned equation is constant, it remains unchanged as the parametric values vary. Therefore, only the first two terms of the equation can be minimized. Hence, the fitting problem can be expressed as: Xðt; PÞ � 0: However, this minimization problem provides the desired fitting with feasible parameter values if P is identifiable, which was confirmed in two steps. First, the structural identifiability was examined, and next, the practical identifiability was confirmed.
Structural identifiability refers to the existence of a unique solution to X(t, P) for each P under initial conditions. If any component of P is implicitly related, different values of P may yield the same solution X(t, P) for a given initial condition, which may hinder the unique estimation of parameter P from the data.
The Fisher information matrix (FIM) is used to confirm the structural identifiability. For a set of observations at n distinct points, a system of 5-dimensional state vector, and 4-dimensional vector of parameters P = (p 1 , p 2 , p 3 , p 4 ) = (β, β f , μ f , I(0)). Thus, the sensitivity matrix S, consists of n time-dependent 1×4 blocks, Aðt i Þ. Here, Here, S is called the sensitivity information matrix (see [36] for more details) and the 4×4 According to the aforementioned definition, the FIM for our problem consists of four columns corresponding to the parameters to be estimated. S is evaluated at P 0 and the small perturbation about P 0 ¼ �0:01P is denoted by DP. HereP denotes the estimated parameter values, that is, DP ¼P À P 0 ¼P � 0:01P. This local perturbation gives rise to a small perturbation, DX ¼ Xðt;PÞ À Xðt; P 0 Þ. Next, the chain rule of differentiation is applied to obtain DX ¼ SDP, indicating that X(t, P) is structurally identifiable if DX ¼ SDP has a unique solution for DP. This is only possible when the FIM rank is equal to the number of unknown parameters [36]. Because the rank of the FIM in all three waves is four, which is equal to the number of unknown parameters, the structural identifiability of the parameters in all three waves is ensured.
Practical identifiability, or ability to be estimated, refers to the sufficiency of available observations, as too few observations may be insufficient for fitting. To evaluate practical identifiability, the profile likelihood of parameters β, β f , and μ f is computed along with the initial value I(0) for three different waves. Profile likelihood revealed the dependency of the NLF on individual parameters, which in turn helped determine the finite confidence intervals for each parameter; otherwise, practical non-identifiability is proved. The related profile likelihoods can be defined as NLFðPÞ; where p i 2 ½p i ð1 À 0:01Þ;p i ð1 þ 0:01Þ�: To determine the confidence interval, we consider 2ðNLFðPÞ À NLFðPÞÞ � w 2 4 . Here, w 2 4 is representing the 95 th percentile range of the chi-square distribution with a degree of freedom equal to the number of unknown parameters. Therefore, the NLF threshold for the 95% confidence interval is NLFðPÞ þ 9:4877=2.
Figs 2-4 represent the maximum likelihood fitting of the time series data for Korea, Pakistan, and Japan, respectively. While the graphs of relative error for Korea, Pakistan, and Japan are given in Fig 5. The estimated parameters are listed in Table 2, and the estimation of the infectious population, I(0), is presented in Tables 3-5 along with the 95% confidence intervals. In addition, the estimated parameter values are biologically meaningful.
Figs 6-8 depict profile likelihoods, which reveal that NLF is minimized at the estimated parameter values, thus confirming practical identifiability.

Results
A compartmental modeling approach is implemented with two different classes of susceptible individuals with different transmission rates. Individuals maintaining strict preventive         measures were grouped as behavior-changed (aware) susceptible, S f . It is assumed that susceptible individuals from S compartment moved to S f compartment upon realizing the high number of active cases or when the preventive measures (as a consequence of the high number of active cases) were strengthened by the public health authority at a rate β f . A behavior-changed individual moves from S f to S following individuals in the S and R compartments at a rate μ f . By analyzing the dynamics of the model, it is confirmed that an outbreak occurs when the basic reproduction number R 0 is above unity. It is worth noting that R 0 increases with transmission rate β (Eq (8)). In addition, higher values of β are associated with a higher and quicker appearance of a peak of the active cases (please refer to Fig 9). The effect of the behavioral change on the intensity and timing of the peaks is discussed in the next subsection. Further, we have fitted our model to time-series data from South Korea, Pakistan, and Japan to understand the evolution of individual responses to public health guidelines.

Behavioral response to public health guidelines modulates the peak
The intensity of the peak is modulated by the behavioral responses of individuals. Figs 10-12 show the peaks as a function of β f and μ f for the three peaks in Korea, Pakistan, and Japan, respectively. The surface plots show that for a fixed β f , the peak increases when people are negligent (i.e., μ f increases) in following the public health guidelines. The intensity of the peak is more sensitive to β f than to μ f . In addition, the sensitivity is higher for the set of values of β f and μ f , lying between the blue and yellow contour lines in the β f −μ f plane. For each value of μ f ,  if the value of β f is above a threshold (shown by the blue contour line), the peak is minimized; that is, maintaining a high behavior change rate is crucial for flattening the next outbreak.

Change in behavioral response during the pandemic: Case study of Korea, Pakistan, and Japan
According to [37], people indicated a minimum of one month and a maximum of two years (mean = 7.20 months) for a return to normal life. The proposed model was used for parameter estimations to examine the effect of individuals' behavioral responses to prolonged public health restrictions established in response to the COVID-19 outbreak. According to our estimation the basic reproduction number R 0 , is 1.4474 (1.5861, 2.1237), 2.5777 (4.4701, 3.3268), and 3.9275 (7.1572, 5.2085) for the first, second, and third waves in Korea (Pakistan, Japan), respectively. Consistently increasing R 0 , in successive waves is a consequence of the appearance of highly transmissible mutants [38,39]. Our estimates of β (Table 2) for the consecutive waves also followed the increasing transmissibility of the new strains. However, the intensity of the peaks did not follow the transmission rates (Figs 2-4). Each country has implemented different strategies to control the disease. Public health authorities implemented necessary actions such as screening, placing people in quarantine after confirmation, closing schools and offices, and making it compulsory to wear a mask in public places when the cases seemed to be high. As a result, in Korea, Pakistan, and Japan, β f increased from the first to the second wave, indicating that more individuals became motivated to strictly follow public health guidelines, outweighing the effect of the increased transmission rate (β) and resulting in a smaller second wave than the first wave. It is worth noting that in Korea, Pakistan, and Japan, μ f increased from the first to the second wave. As the intensity of the peak is more sensitive to β f than to μ f (as mentioned in the previous subsection), the effect of β f dominates. Additionally, during the third peak, β f decreased in Korea and Japan but remained constant in Pakistan (comparing columns 4 and 5 of Table 2). The third peak is the greatest in Korea and Japan because of poor awareness and strong transmissibility. Although the third peak in Pakistan is larger than the second in the country because of the increased transmissibility, constant awareness maintained it below the first peak.
For the estimated value of transmission rate in each wave, we plot the peaks of the active cases with respect to β f and μ f in Figs 10-12. These figures show that the peak of the active cases is lower when the values of β f are high, indicating that when fear of the disease is high, the peak can be delayed and reduced. In the case of a higher behavioral change ease rate, the peaks are higher, and vice versa. This indicates that a reduction in the fear of the disease may cause more transmission in society and may cause the next wave. If the transmission rate is high, then, high fear of the disease would reduce the number of active cases.

Discussion and conclusion
In this study, a mathematical model is used to explore the role of public acceptance of the public health guidelines during a pandemic outbreak. It is challenging to maintain a high level of awareness of preventive measures during a long-term disease with the emergence of subsequent outbreaks [37]. Therefore, a model that differentiates between susceptible and aware susceptible individuals in terms of disease transmission is proposed. Our analysis revealed that apart from the transmission rate and basic reproduction number, the behavior change rate modulates the intensity and timing of the peak. It is assumed that when individuals know their infection status, they are either hospitalized or isolated at home. Therefore, the infected class is not divided based on behavioral changes, as in the case of susceptible individuals. The fear of the disease causes susceptible individuals to move to the behavior-changed susceptible class, or because of monotonous maintenance of preventive measures, the behavior-changed susceptible individuals move back to the susceptible class. Our mathematical analysis confirmed that although R 0 is the pivotal threshold for the COVID-19 outbreak, behavioral changes associated with preventive measures may have a considerable effect on the course of the pandemic.
Typically, data from active and recovered cases are used to fit the model, and the infectious compartment handles the spread of the disease. However, reliable data for the infectious class are not currently available because of a significant number of asymptomatic infectious individuals. To overcome this limitation, the initial values of the infectious class were estimated using the maximum-likelihood method. For a comprehensive analysis of the model, an identifiability analysis is performed to guarantee that the estimated unknown parameters were uniquely determined and that the numerical values of the estimated parameters were therefore meaningful.
Different techniques have been used to fit data with mathematical models, several of which have been described in [40][41][42][43]. In this study, the maximum likelihood method is used for fitting. As case studies, the COVID-19 outbreak in Korea, Pakistan, and Japan was investigated to determine how individual tendencies to follow preventive guidelines may affect the outbreak. The transmission, behavior change, and behavioral change ease rates were estimated for the three waves of the COVID-19 outbreak in South Korea, Pakistan, and Japan. The behavioral change rate reached an all-time high at the time of the second wave in Korea and Japan owing to government-imposed restrictions, such as school closures, public mask-wearing, and social isolation. However, it is not maintained in the third wave; therefore, the peak in this wave is the highest among the active cases.
Pakistan lacked diagnostic capabilities during the early stages of COVID-19; thus, suspicious samples were submitted to overseas labs [44]. Even though the transmission rate was the lowest during the first wave despite the low behavioral change rate and lack of testing facilities, the peak of this wave is the highest. Pakistan is then given test kits from China as well as primers from Japan [44]. Pakistan could examine samples from suspected cases nationwide. During the second wave, the fear level increased substantially. Compared to the second wave, the transmission rate of the third wave increased by 70%. However, the fear level is approximately the same as in the second wave. Therefore, Pakistan's third peak is higher than its second.
It should be noted that a cross-country comparison is not carried out, as the country-based size of the peak and behavioral response would vary depending on the respective population size and socioeconomic conditions [45]. Rather, different peak scenarios were compared for each country. It is observed that in Korea Pakistan and Japan, individual responses to public health guidelines evolved with disease progression, and people resumed normal life, perceiving the risk as having declined. Usually, the fear of the disease prompts people to respond by limiting their contact; however, when they face insufficient incentives to alter their behavior, they return to normal behavior, which may modulate the size of the upcoming outbreak. Our mathematical findings highlight the critical necessity of behavioral change educational interventions, and community-level media campaigns. To address these issues, policymakers may improvise their strategy to make people more inclined to follow safety instructions [46].